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Abstract We introduce a three-dimensional, computationally feasible, mesoscopic model for 
^ snow crystal growth, based on diffusion of vapor, anisotropic attachment, and a semi-liquid 

^ boundary layer. Several case studies are presented that faithfully emulate a wide variety of 

CS| physical snowflakes. 
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1 Introduction 



In this paper we exhibit some virtual snowflakes, or snowfakes, generated by a natural, fully 
three-dimensional algorithm for snow crystal evolution. The present study extends our earlier 
work on growth and deposition |GGH [GG2^ IGG3] . and other previous efforts in this direction 
|Pac[ IReij . The key features of our model are diffusion of vapor, anisotropic attachment of 
water molecules, and a narrow semi-liquid layer at the boundary. All three ingredients seems 
to be essential for faithful emulation of the morphology observed in nature. The algorithm 
assumes a mesoscopic (micron) scale of basic units for the ice crystal and water vapor, which 
eliminates inherent randomness in the diffusion and the attachment mechanism. This brings 
the process within reach of realistic simulation; by contrast, any three-dimensional approach 
based on microscopic dynamics is completely beyond the scope of present computing technology. 
We refer the reader to |GG3j for a brief history of snow crystal observation and modeling, 
background on our approach in a two-dimensional setting, and many references to the literature. 
See also [NRj for another attempt at spatial mesoscopic modeling. 

There are many papers and books, for a variety of audiences, dealing with snowflake pho- 
tography and classification, the underlying physics, or some combination thereof, so we will not 
offer a comprehensive review here. Excellent introductions to the subject include the classic 
book by Nakaya [Nakj . early empirical studies and classification schemes |BH| and |ML| . and 
more recent papers and books by K. Libbrecht [Libl] ILib2[ ILib3[ ILib4[ ILib5[ ILRj . Among 
research papers that attempt to decipher the three-dimensional aspects of snow crystals, the 
standout reference is [TEWFj; also worth mentioning are |Iwaj . |NKj and |Nel| . The single 
most convenient resource for comparison of our simulations to physical crystals is Libbrecht 's 
field guide [Lib6]. 




Fig. 1. Tip instabihty and oblique top (left) and bottom (right) views of the final crystal. 
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As a preview of the capabilities of our model, let us illustrate the crystal tip instability and 
initiation of side branching studied in the laboratory by Gonda and Nakahara [GNj . A sequence 
of four still frames from their paper was reproduced in [GG3] so we will not show it here. But 
Fig. 1 depicts the top view of a corresponding snowfake at four different times (12, 15, 18, and 21 
thousand), and oblique views of the crystal's top and bottom at the final time. The parameters 
are: /3oi 2.8, (3io f320 2.2, (3u /32i 1.6, /^ao = /^ai = 1, = -005, /xio M20 -001, 
fi = .0001 otherwise, = .01, and p = .12. Their role, and that of the initial state, will 
be described in Section 2. Similarity between the real and simulated sequences is striking: in 
both instances a defect arises at a characteristic distance from the crystal tip, becomes more 
pronounced, and later gives rise to a side branch with its own ridge structure similar to that of 
the main branch. Note also that our snowfake has its ridges and most of its markings on the top 
side; the bottom is almost featureless. This is due to a small downward drift in our model, an 
aspect we will discuss later in more detail. The direction of the drift represents the motion of the 
crystal in the opposite direction — we prefer upward motion because interesting features then 
appear on top, although this would obviously correspond to the bottom of a falling snowflake. 
We should also note that the drift value means that, during its evolution, our simulated crystal 
moved for about 200 space units, which is comparable to the diameter it reached. This is 
typical of drift values that erase features on one side without otherwise significantly changing 
the morphology. Our model thus predicts that a significantly larger range of motion during 
growth is not possible for most interesting physical snow crystals, such as dendrites or plates. 
Another example of our algorithm's potential to make new predictions about basic aspects of 
snow crystal growth is the location of markings. From micrographs, it is almost impossible to 
tell whether these are on the top, bottom, or inside a given physical specimen, so little attention 
has been paid to this issue to date. We have gathered a considerable amount of evidence that 
inside markings are quite common (cf. Sections 7, 8 and 9). 

Our account will focus on seven case studies that reproduce many features commonly ob- 
served in actual snowflakes: ridges, ribs, flumes and other "hieroglyphs," formation of side 
branches, emergence of sandwich plates, hollow columns, hollow prism facets, and so forth. We 
also explore dependence on the density of vapor, and the aforementioned effect of drift, and inhi- 
bition of side branches by the semi-liquid layer. Varying meteorological conditions during growth 
are considered very important |Lib6| so we include several examples, such as plates with den- 
dritic tips and capped columns, that are believed to arise due to sudden changes in the weather. 
However, we will encounter snowfakes that grew in a homogeneous environment but give the 
impression that they did not. We will occasionally address dependence of the final crystal on its 
early development, and conclude with a few eccentric examples that may be too brittle to occur 
in nature. These typically arise near a phase boundary, when the dominant direction of growth 
is precarious. A complete collection of snowfakes from our case studies (with some additional 
information, such as simulation array sizes), and a slide show are available for download from: 

http : //psoup . math . wise . edu/Snowf akes . htm 

The first order of business, in the next section, is to describe the snowfake algorithm in detail. 
Four subsequent sections discuss computer implementation and visualization tools, mathematical 
foundations, parameter tuning, and extensions of the model. The remainder of the paper is then 
devoted to the case studies. 
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2 The algorithm for three-dimensional snow crystal growth 

Our basic assumptions cirG clS follows: 

Al. The mesoscopic (micron-scale) building blocks are (appropriately scaled) translates of the 
fundamental prism, which has hexagonal base of side length 1/a/3 and height 1; 

A2. In its early stages of growth, from microscopic to mesoscopic, the crystal forms a hexagonal 
prism, and then it maintains this simple polyhedral shape until it reaches the size of a few 
microns across. 

A3. Diffusion outside the growing crystal is isotropic except possibly for a small drift in the 



A4. Crystallization and attachment rates depend on the direction and local convexity at the 
boundary; 

A5. There is a melting rate at the boundary, creating a quasi-liquid layer. 

Note that the side (rectangular) faces of the fundamental prism are commonly referred to as 
prism faces, while the top and bottom (hexagonal) ones are called basal faces. 

The lattice for our model is T x Z, where T is the planar triangular lattice (see Fig. 2). 
This is not precisely the crystalline lattice of hexagonal ice Ih, which is obtained by removing 
certain edges and sites from T x Z, and then applying a periodic deformation [NRj, but we are 
constructing a mesoscopic model that should obscure such fine details. Therefore, each x G T x Z 
has 8 neighbors, 6 in the T-direction and 2 in the Z-direction. 

At each discrete time t = 0, 1, 2, . . . and with each site x G T x Z, we associate a Boolean 
variable and two varieties of mass: the state of the system at time t at site x is ^t{x) = 
{at{x),bt{x),dt(x)) where the attachment flag 



bt{x) = the boundary mass at x at time t (frozen if at{x) = 1, quasi-liquid if at{x) = 0), 



Our dynamics assumes that the diffusive and the quasi-liquid mass both change to ice when the 
site joins the crystal, and stay in that state thereafter. The two types of mass can coexist on 
the boundary of the snowfake, but only boundary mass persists inside the snowfake while only 
diffusive mass occurs outside and away from the boundary. 

The initial state will consist of frozen mass 1 at each site of some finite set, on which also 
ao = 1, with ao and bo = and do = p everywhere else. In keeping with assumption (A2), the 
most natural choice for this finite set, a singleton at the origin, often does not work well, as its 
Z-direction neighbors see 7 neighbors off the crystal's boundary. This means that it is common, 
even for low p, that the dynamics immediately triggers a rapid expansion in the Z-direction. To 



Z-direction; 




if X belongs to the crystal at time t 
otherwise; 



and 



dt(x) = the diffusive mass at x at time t 



(vapor). 
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prevent this singularity, our canonical initial state consists of a hexagon of radius 2 and thickness 
1, consisting of 20 sites. Other non-symmetric initial states will be discussed later. 

Let us now describe the update rule of our snowflake simulator, which performs steps {i)-{w) 
below in order every discrete time unit. The reader should observe that total mass is conserved 
by each step, and hence by the dynamics as a whole. 

Write A/J = {x} U {y : ?/ is a neighbor of x in the T-direction}, — {x} U is a 

neighbor of x in the Z-direction} for the T- neighborhood and Z- neighborhood of x, respectively. 
We also let J\fx=Mj UAf^, and set 

At = {x : at{x) = 1} = the snowfake at time t; 

dAt = {x ^ At : at{y) = 1 for some y G Mx} = the boundary of the snowfake at time t] 
At = AtU dAt. 

The complement of a set A is denoted by A^. Also, we use ° (degree) and ' (prime) notation to 
denote amounts of mass before and after a step or substep is completed. If there is more than 
one intermediate step, we use double primes. This is necessary since some mass allocations may 
change more than once during a single cycle of the steps. At the end of each cycle the time t 
advances to t + 1 . 

Steps of the update rule: 

i. Diffusion 

Diffusive mass evolves on A^ in two, or possibly three, substeps. The first substep is by 
discrete diffusion with uniform weight ^ on the center site and each of its T-neighbors. Reflecting 
boundary conditions are used at the edge of the crystal. In other words, for x G A^, 

(la) <{x) = \ Yl dtiy)- 

The second substep does the same in the Z-direction: 

(lb) <'(x) = ^4(x) + A <{y)- 

For X G dAt any term in the sum in (la) (resp. (lb)) corresponding to y ^ At is replaced by 
d^{x) (resp, d[{x)). 

The reason for the weights in (lb) is as follows. Imagine we tessellate with translates of 
the fundamental prism and scale the lattice T x Z so that the lattice points are in the centers 
of these prisms. The "bonds" in the top left frame of Fig. 2 thus all have unit length and we 
eventually visualize the crystal by drawing prisms that are centered about sites of At. Rule (lb) 
ensure that diffusion on the scaled lattice is isotropic, in agreement with assumption A2. 

As mentioned in the Introduction, there is also good reason to consider the more general 
case of diffusion with drift in the Z-direction, corresponding to downward (or upward) motion 
of the snowflake. The third diffusion substep is thus: 

(Ic) d'l\x) = (1 - . (1 - at{x - es)) • d'l{x) + • (1 - at{x + 63)) • d'l{x + 63), 
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where 63 = (0,0, 1) is the third basis vector. Parameter (j) measures the strength of the drift, 
and needs to be smah for the dynamics to remain diffusion-hmited. 

ii. Freezing 

Assume that x G dAf, and denote 

(2a) nj{x) = #{y G A/^ : a°{y) = 1} A 3, nf{x) = #{y G Af^ : a°(y) = 1} A 1. 

Proportion 1 — nf (x)) of the diffusive mass at x becomes boundary mass. That is, 

h[{x) = + (1 - K{nJ{x), nf{x)))d°{x), 
d[ix) = KinJix),nfixMix). 

The seven parameters i G {0, 1}, j G {0, 1, 2, 3}, i+j > 0, constitute one of the ingredients 

that emulate the dynamics of the quasi-hquid layer at the boundary of the crystal. The other 
ingredient, /x, appears in step iv below. We assume that decreases in each coordinate since 
"more concave corners" at the boundary dAf, i.e., those with more neighbors in At, should catch 
diffusing particles more easily. 

Hi. Attachment 

Assume again that x G dAt and define the neighborhood counts as in (2a). Then x needs 
boundary mass at least f3{nj {x),nf{x)) to join the crystal: 

(3) If bt{x) > (3{nJ{x), nf (x)), then at{x) = 1. 

Again, we have seven parameters f3(i,j), i G {0, 1}, j G {0, 1, 2, 3}, i+j > 0, and the assignment 
only makes physical sense if (3 decreases in each coordinate. 

In addition, we assume that a^(x) = 1 automatically whenever nj{x) > 4 and nf{x) > 1. 
This last rule fills holes and makes the surface of the crystal smoother, without altering essential 
features of the dynamics. 

At sites X for which aj(x) = 1, the diffusive mass becomes boundary mass: bt{x) = 6^(x) + 
(i^(x), dt{x) = 0. Attachment is permanent, and there are no further dynamics at attached sites. 
Thus we do not model sublimation, although it may play a significant role in the last stages of 
snow crystal evolution (cf. p. 27 of |Lib6| ). 

iv. Melting 

Proportion /x(nj^(x), nf (x)) of the boundary mass at each boundary site becomes diffusive 
mass. Thus, for x G dAf, 

b',{x) ^ {1 - f,{nj{x),nf{x)))b°{x), 
d'tix) = d°ix) + ^,inJ{x),nf{x))b°{x). 

Again, /i is decreasing in each coordinate. 

Fig. 2 summarizes our model in three frames. At the upper left is a portion of the underlying 
lattice T X Z. The central site represented as a larger black ball has its neighborhood indicated 
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in black, and a translate of the fundamental prism is centered at that site. In the upper right 
detail, blue translates of the fundamental prism are drawn around each site of a small crystal. 
Seven boundary sites are depicted in red and each is labeled by its boundary configuration. 
For example, the "21" site has 2 horizontal (T-) neighbors and 1 vertical (Z-) neighbor, and 
consequently needs boundary mass /32i to join the crystal. Finally, the lower panel shows a 
flowchart for the algorithm. There are three epochs in the life of a site. Away from the crystal's 
boundary, it only exchanges diffusive mass dt with its neighbors. Once the crystal grows to reach 
the site's neighborhood, two additional effects, melting and freezing, promote exchange between 
diffusive mass dt and boundary mass bt. Final changes occur once boundary mass exceeds the 
threshold (3 (which depends on the neighborhood configuration): the site attaches and the two 
types of mass merge into bt. 





non-boundary 




boundary 








dt 


k 


K, 








bt 





attached 




Fig. 2. The stacked triangular lattice T x Z {top left), coding of boundary configurations 
{top right), and a flowchart for the growth algorithm {bottom). 



3 Notes on computation and visualization 

Following the same strategy as for our previous two-dimensional model |GG3| , the dynamics 
actually run on the cubic lattice Z^, which can be mapped onto x Z. Our basic computa- 
tional engine is written in C, but MATLAB is used for mapping and visualization. As mentioned 
previously, the snowfakes are depicted by drawing visible boundaries of translates of the funda- 
mental prism centered on sites of At. Since this straightforward procedure makes jagged vertical 
boundaries, we apply a smoothing algorithm at the boundary that enlarges the crystal by no 
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more than one mesoscopic unit. (This algorithm is not apphed to the smah snowfake in Fig. 2.) 
MATLAB's patch routine renders the faces. For better results we then emphasize edges using 
the line routine. 

MATLAB's visualization tools certainly provide adequate representations for detailed inves- 
tigation of the resulting crystals. They do not, however, give a satisfactory comparison with the 
best snowflake photographs |LR^ ILibSj ILibGj . typically taken from directly above the (predom- 
inantly two-dimensional) crystal, which is in turn illuminated from below. This viewpoint can 
be effectively simulated by ray-tracing, as implemented here by the POV-Ray software jPOVj . 
Our program automatically outputs a file with a triangulation of the crystal's boundary, which 
is then used by the mesh2 command in POV-Ray. 

We would like to point out that both the algorithm and visualization procedures require 
considerable computing power and memory. At present (fall 2007) , our simulations are very time 
consuming, barely feasible on commercial personal computers. (In fact, an adaptive resolution 
algorithm is necessary to make the boundary descriptions manageable.) Progress in studying 
snowfakes is therefore quite slow, precluding systematic classification of the dynamics. Our 
goal has been to find representative examples that seem to replicate physical snow crystals and 
thereby shed light on their evolution. 

For computational efficiency, if the diffusion step is isotropic one can exploit symmetry by 
taking the finite lattice to be a discrete hexagonal prism with patched wrap edge conditions. 
When = and the initial state has complete symmetry, it thereby suffices to compute the 
dynamics on ^ of the whole space. There are two good reasons for giving up complete symmetry 
of the rule. First, the initial state may not be symmetric, and second, the diffusion may have a 
drift. For computational efficiency, we only give up reflectional symmetry around the x?/-plane 
(recall that the drift is only in the Z-direction), allowing the initial state to depend on the z 
coordinate, but retaining its hexagonal symmetry in the x and y coordinates. This increases the 
space and time demands of the fully symmetric program by a factor of 2. 

The program stops automatically when the density at the edge of the lattice falls below a 
given proportion of the initial density (typically 2p/3 or p/2), or when the crystal gets too close 
to the edge (snowfake radius greater than 80% the radius of the system). 

4 Connection to pde, and size of the parameter space 

Mathematically, our algorithm is a discrete space and time version of a free boundary^ or Stefan^ 
problem |Lib2[ ILib3[ ILib4| . This is a partial differential equation (pde) in which the crystal 
is represented by a growing set At and the density (i.e., supersaturation) of vapor outside it as 
u = u(x^t). Then i/, is on the boundary dAt, and satisfies the diffusion equation outside the 
crystal 



The velocity of the boundary at a point x G dAf with outside normal u is given by a function 



(1.1) 



— = An, xe Al 



(1.2) 
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Considering the slow growth of At^ diffusion equation (1.1) may be simphfied to its equihbrium 
counterpart Au = |Lib2[ ILib3[ ILib4| , which makes this into an anisotropic version of the 
Hele-Shaw problem. 

Presumably under diffusion scaling, in which space is scaled by e, time by e~^, and e ^ 0, 
the density field and the occupied set in our model converge to a solution of the Stefan problem. 
We hope to provide rigorous justification for this connection, and identification of the limit w 
in terms of model parameters, in future work. 

The boundary velocity function w = w{X, v) is defined for A > and three-dimensional unit 
vectors u G S^. In order to develop a rigorous mathematical theory, the most convenient assump- 
tions are that w is continuous in both variables, nondecreasing in A, and satisfies w{X^ v) < CX 
for some constant C independent of A and u. Under these conditions, the non-isotropic Stefan 
problem (1.1-1.2) has a unique viscosity solution at all times t > 0, starting from any smooth 
initial crystal. This is proved in [Kimj for the isotropic case (when w is constant); assuming 
the listed properties of the proof extends to our general setting. We should note, however, 
it has long been known that the crystal's boundary will not remain smooth [SB]. Indeed, this 
will be no mystery once we present our simulations, which feature a considerable variety of 
singularities and instabilities. Presumably these make direct numerical computation with the 
pde very challenging, explaining why numerical pde-based models for snow crystal growth have 
not been satisfactory (cf. |Sch] ). For further mathematical theory and references, we refer the 
reader to jKim^ ICK| . 

For the model studied here, w{X,u) will be linear in A, since the attachment and melting 
rates are independent of the vapor density. This may not always be the case; in fact, some of 
the literature even considers the possibility that w is non-monotone in A [Lib3l IGG3j . Analysis 
of such cases would present new theoretical challenges, and from simulations of our 3d model 
it appears that nonmonotonicity is not needed for observed phenomena in nature. Monotone 
nonlinearity, arising from monotone density dependent rates, is harder to dismiss and worth 
further investigation - for instance, it is possible that w vanishes for very small A. 

Once we accept that our scheme approximates the viscosity solution of (1.1-1.2), the macro- 
scopic evolution of the crystal is uniquely determined by its initial state and the velocity function 
w. In turn, w is determined by very few physical parameters, perhaps just two: temperature 
and atmospheric pressure jLib2[ ILib3^ ILib4]. Therefore, possible evolutions from a fixed seed 
comprise a three-dimensional manifold (its coordinates being the supersaturation level, tem- 
perature, and pressure) in an infinite-dimensional space of possible velocities w. Much of the 
ongoing snow crystal research constitutes an attempt to understand the structure of this mani- 
fold, a daunting task since the underlying (perhaps quantum) attachment physics is very poorly 
understood, controlled homogeneous environments are hard to design, and crystal evolution is 
difficult to record. Our model does not have these problems. Instead, its main weakness is the 
number of free parameters that need to be tuned to approximate it; at a particular temperature 
and pressure. It helps that our parameters have intuitive meaning, but finding a particular 
realistic snowfake involves approximating an a priori infinite-dimensional object w by one of 
finite but high dimensionality. The challenge is compounded by very incomplete information - 
all that is typically observable in nature is the final crystal, which may have been subjected to 
numerous changes in conditions and orientation during growth, as well as sublimation and per- 
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haps even artifacts of the recording process. It is thus no surprise that our parameter selection 
is an arduous and imprecise task. 

In the next section we wih describe some ad hoc rules that we have used to generate our case 
studies, but the issue of parameter selection is in dire need of further investigation. What we can 
say is that the best examples are quite sensitive to perturbations in w. Thus they require good 
approximations and a large number of judicious parameter choices. In addition, the dependence 
on the initial seed is often quite dramatic. These observations underscore both the marvel and 
the fragility of natural snowflakes. 

At the same time, we wish to emphasize the conceptual simplicity of our model. The large 
parameter space is a consequence of geometry rather than an excessive number of modeling 
ingredients. Apart from the two scalar parameters - density p and drift - we have only three 
vector parameters — attachment threshold freezing rate 1 — and melting rate p — whose 
high dimensionality arises from the many possible boundary arrangements. The parameter set 
can be reduced, but some tuning will always be necessary, as illustrated by the "random" crystal 
in Fig. 3. This was obtained by choosing Hi = .1, fi = .001, p = .1, = 0, and all /3's equal 
to 1 except (3oi = 1.73 and /3io = (^20 = 1-34. These values are in a sensible neighborhood of 
the parameter space, but the last two attachment rates were selected by chance. The result 
has some physically reasonable features, but one immediately notices an excessive density of 
branches and inordinately high ridges. 



Fig. 3. A "failed" snowfake. 



5 Effective choice of parameters for simulations 

While optimal choices of parameters requires considerable guesswork, there are a few guidelines 
we have developed. Some come from mathematical arguments, others from experimentation; 
both are described in this section. 

Our simulator represents diffusion by discrete averaging in time t, which is also discrete. The 
bulk effects of this operation expand at the rate y/i, although the extreme radius of its influence 
(or light cone) grows linearly in t. If the initial density p of our discrete vapor field is too large, 
then the crystal may expand in some direction as fast as the light cone, or perhaps fall behind it 
by 0{y/i). We call parameter sets leading to this behavior the Packard regime; it is clearly not 
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physical, as it depends on the discrete nature of the averaging. However, systems of this sort are 
able to generate fractal plates reminiscent of Packard snowflakes [PacJ IGG3 ] and exhibit one 
variety of faceting (cf. [NR]). In our simulations we systematically avoid the Packard regime by 
keeping the density low. For the extremal points of our snowfakes not to expand at light speed, 
the conditions are 

(1 - A^oi)p < /?01, (1 - I^10)p < Ao, 

as is easy to see from the description of the rule. Our densities are typically considerably smaller, 
since large densities generate expansion that is too rapid to be realistic, at least in its initial 
stages. As mentioned previously, a surprisingly important role is also played by the choice of 
initial seed. 

On the other hand, it is clear that a very large melting rate will stop growth altogether. This 
happens if the flow out of the boundary mass exceeds the flow in just before that mass exceeds 
the threshold for attachment. A sufficient condition for continual growth in all directions is 
therefore 

Moi/5oi < (1 - ^oi)p, MioAo < (1 - /^io)/>, 

since the 01 and 10 boundary arrangements always have the slowest potential growth. In the 
great majority of examples we will present, parameters for the 20 and 10 arrangements agree. 
In this case, the last condition is necessary as well — if it does not hold, then the growth is 
convex-confined in the T-direction. 

Let us now describe a few rules of thumb when searching for snowfakes that emulate nature. 
We commonly start with a reduced parameter set. Namely, we set the k^s to a common value, 
say, = .1. Then we select two different /3 parameters, /3oi and /3io = /320 — f^iii with all the 
remaining /3's fixed to 1. The size of P20 controls the strength of the convexifying mechanism, 
assumed to be the same in both the xy and z directions. Indeed, if /32o is large, then the 
crystal will remain a perfect hexagonal prism for a long time. The only other parameters are 
the common value of all /x's and the vapor density p. This is a more manageable four-parameter 
space that encodes four essential elements of three-dimensional snowflake growth, each with a 
single tunable parameter: diffusing supersaturation level (p), convexifying strength {[^20), quasi- 
liquid layer smoothing (/x), and preference for the Z-direction over the T-direction (/3oi//32o)- 
This scheme is used to identify the neighborhood of a desired morphological type in phase space. 
Then parameters are perturbed for added realism. 

One of the most important lessons of our two-dimensional model |GG3j was that the melting 
parameter p inhibits side-branching and is therefore important for dendrite formation. When 
/i = 0, it seems impossible to avoid an excessive density of branches. Indeed, this role of /i is 
easily understood. Namely, fi creates a positive density at the boundary, due to flow out of 
the boundary layer. This density has the effect of reducing the ambient vapor density by a 
fixed amount, independent of location, and hence disproportionately affects regions of smaller 
density. (To a very rough first approximation pLib4], the expansion speed is proportional to 
■sfpl\ft when /i = 0.) Since there is clearly less mass between branches than at the tips, growth 
and side branching there gets stunted by increasing [x. 

Realistic "classic" dendrites occur for a relatively narrow range of choices for /x, once the 
other parameters are held fixed. Typically, though, the other parameters need to be perturbed 
along with /x; increasing [x alone tends to erode all complex structure. 
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The markings seen on snow crystal plates are sometimes called hieroglyphs. These often 
have fairly regular geometric forms, such as ridges, flumes, ribs, and circular shapes, but can 
also exhibit more chaotic patterns. In photomicrograph collections |BH[ ILR[ ILib5|. ILib6| it 
is usually unclear whether the marks are on the outside of the crystal or within what we call 
sandwich plates. In our experiments, the inner structures are much more prevalent, so we are 
glad to observe that they are abundant in nature [EMPJ. To obtain nice outer markings, the 
ratio /3oi//32o needs to be sufficiently large, but there is then a tendency for the crystal to become 
too three-dimensional. Again, the correct choice is often rather delicate. Inner markings occur 
generically for small values of this ratio. 

Finally, different /^'s may appear to be a more natural mechanism to enforce anisotropy 
than different /3's, as they directly correspond to sticking, or killing, of particles at the crystal's 
boundary. However, for this effect to be significant, the /^'s need to be very close to 1; otherwise 
the killing at the crystal boundary is too rapid to make a difference, and then the already slow 
growth proceeds at an even more sluggish pace. While less physically appealing, we view the 
/3's as a reasonable compromise for the sake of computational efficiency. 

6 Variants and extensions of the model 
6.1 Uniform snowfakes 

Since attachment thresholds /3 vary, the mass of the final crystal is not uniform. There is 
a variant of our algorithm that removes this defect with little change in observed morphology. 
Assume that there is no automatic ffiling of holes; instead, boundary mass exactly 1 is needed for 
attachment when nj{x) > 4 and nf{x) > 1. Then a uniform crystal is obtained by performing 
the following additional step just after step Hi in the simulator: 

iii\ Post-attachment mass redistribution 

To redistribute any excess mass from the attached site to its unattached neighbors, let 



6.2 Simulation without symmetry 

As explained in Section 3, at the cost of a 24-fold slowdown compared to our fully symmetric 
model, implementation of the algorithm without exploiting symmetry makes it possible to study 
the evolution from arbitrary initial seeds. Such an extension is necessary in order to produce 
snowfakes corresponding to exotic forms such as triangular crystals, split stars, and bullets. We 
have conducted a few experiments along these lines with our planar model [GG3J . but in three 
dimensions a simulator dramatically faster than our current one is needed. We have future plans 
to develop a suitably high-performance parallel version. 



be the number of non-attached neighbors. Then, for every x with a^{x) = 0, 
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6.3 Random dynamics 

Our only three-dimensional snowfakes to date are deterministic, since randomness would also 
require the just discussed simulation without symmetry. We propose to include an additional 
parameter e representing residual noise on the mesoscopic scale, as we did in the two-dimensional 
setting [GG3]. Again, e would need to be quite small, say on the order 10~^. The random 
perturbation of diffusive mass from [GG3j is not suitable in 3d since it is not physical to violate 
mass conservation. Instead, a small random slowdown in the diffusion rate is more appropriate. 
To this end, first denote the (hnear) operation on the field in (la-lc) by V; thus step i can 
be written as d^^^ = V{d^). Next, let ^t(x), t > 0, x G T x Z, be independent random variables, 
equal to e > or 0, each with probability 1/2. Here the field ^ represents the proportion of 
particles that refuse to diffuse at position x and time t. The randomized step i now reads 



In a natural way, this represents small random temperature fluctuations in space and time. 

Similarly, one could introduce a small proportion of particles that refuse to freeze in (2b), or 
melt in (4); e.g., (2b) would be replaced by 



7 Case study i : ridges and plates 

Our prototypical snowfake has p = .1 and the canonical initial state of radius 2 and thickness 1. 
Fig. 4 depicts the crystal after 70000 time steps, when its radius is about 350. Its parameters 
are /3oi = 2.5, /3io = /520 = f^ii = 2, f3so = /32i = /Jai = 1, = -1, = .001, and = 0. 



< = P((l - 6)d?) + = T^id't) + - md't)' 



h',{x) = 6?(x) + (1 - ^{nj{x), nf{x)M{x){l - ^t{x)). 
4{x) = ,.(n?(x),nf(x))d°(x)(l -6(x)) + d?(x)6(x). 




Fig. 4. The oblique (MATLAB-rendered) and top (ray-traced) views of the crystal. 



7 CASE STUDY I : RIDGES AND PLATES 



13 



We invite the reader to compare the simulated crystal with some of the photographs at 
pLiibS] and especially with Fig. 1(h) in [TEWF], a snowflake obtained at temperature about 
— 13°C. We think of our length unit as about 1/im, so even the sizes of the two objects roughly 
match. Perhaps the most striking features shared by the snowfake in Fig. 4 and physical ones 
are the ridges, elevations in the middle of each main branch, with less pronounced counterparts 
on the side branches. We begin by illustrating how these ridges are formed and maintained. 
In the process we also encounter the branching instability, when the initial growth of a thin 
hexagonal plate is no longer viable and it gives birth to the six main branches. 

As shown in Fig. 5, ridges are formed quite early in the evolution, by mesoscopic bumps 
known as macrosteps that are near, but not too near, the center of the plate. This is how the 
ridges grow (very slowly) in the vertical direction — compare with times 4044 and 7099, which 
also feature such bumps. The ridges spread to a characteristic width, but sharpen to a point 
near the branch tip. One can also observe the commonly observed flumes (called grooves in 
|Lib6| ) that form on both sides of a ridge. 




Fig. 5. The crystal at times 820, 863, 1600, 4044, 5500, 7099, and 9500. 



The small indentation that emerges, due to lower vapor density, in the middle of each prism 
facet at time 5500, has appeared several times before. However, this is the first instance when 
the growth is unable to repair it. Instead, the growth there virtually stops, while the six main 
arms continue to grow and eventually produce two types of side branches: common, relatively 
thick double-plated branches that we call sandwich plates, and more unusual thin plates with 
their own ridges. The tip of each arm assumes its characteristic shape by the final frame of 
Fig. 5. 

It is perhaps surprising how dramatically this scenario depends on the initial (micron scale) 
state. Keeping everything else the same, we change the initial prism to one with radius 2 and 
thickness 3. The previous rather complex and aesthetically pleasing evolution is replaced by a 
growing double plate (Fig. 6). (Remarkably, even adding a small drift does not help matters 
much.) This dichotomy arises frequently in our model — within a neighborhood of the parameter 
space that produces planar crystals there are two stable attractors: one with outside ridges and 
the other a split plate with ridges on the inside. As much of the literature points out, split plates 
are extremely common in physical crystals (cf. [IwaJ ). 
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Fig. 6. Oblique and side views of the crystal from a different initial state. 

Finally, let us experiment with changing the density p. We exhibit five crystals, each with 
the canonical initial condition and all other parameters of the prototype unchanged, but at 
different densities and different final times. Dramatically lower density does promote faceting 
([L ib6[ iLRj ). but a moderate perturbation seems to mainly promote slower growth, without a 
change in morphology. 




Fig. 7. At density p = .15, the side branches have particularly well-defined ridges. 



CASE STUDY I : RIDGES AND PLATES 




Fig. 9. Density p = 



.05 results in sectored plates. 
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Fig. 10. Density p = .045 results in sectored branches. 




Fig. 11. Density p = A results in sandwich plates with inner ridges. 

The example in Fig. 11 (pictured at time 120000) never undergoes the branching instability 
illustrated in Fig. 5, although it does develop fairly standard ridges that persist until about time 
40000. This is the time shown in the first frame of Fig. 12; subsequent frames show the evolution 
in time increments of 10000. We observe that a completely different sandwich instability takes 
place: first the tips and then the sides of the snowfake thicken and develop sandwich plates. It is 
also clear from the time sequence that this morphological change is accompanied by a significant 
slowdown in growth. We should emphasize that this slowdown is not due to the depletion of mass 
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on a finite system: much larger systems give rise to the same sandwich instability well before 
the edge density diminishes significantly. Neither is this slowdown accompanied by a significant 
growth in the Z-direction — in the period depicted, the radius in the Z-direction increases from 
6 to 7, whereas the radius in the T-direction increases from 67 to 87. Instead, much of the 
growth fills space between the ridges, the remnants of which end up almost completely below 
the surface. 




Fig. 12. The crystal of Fig. 11 at earlier times. 



Note that the snowfake of Fig. 10 is also experiencing the sandwich instability at about 
the capture time. The difference in that case is that the growing crystal also experienced the 
branching instability earlier in its development. 



8 Case study ii : classic dendrites 




Fig. 13. p = .105 : a fern dendrite. 

For this series of snowfakes, /3oi = 1.6, /3io = (^20 = 1-5, /^n = 1.4, f3so = (32i = f^si = 1, 
= .1, all /X = .008, = 0, and growth starts from the canonical initial state. We will 
again look at how morphology is affected by different vapor densities p. The simulations argue 
persuasively that the frequency of side branches decreases with decreasing p. When p = .105, 
the branches are so dense that the crystal is rightly called a fern, while the examples with p = .1 
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and p = .095 have the classic look of winter iconography. These are our largest crystals, with 
radii around 400. A more substantial decrease in p eliminates any significant side branching on 
this scale, resulting in a simple star for p = .09. As should be expected from Section 7, further 
decrease finally produces a sandwich instability at the tips, resulting in thick double plates. 
In this instance, slow growth at the branch tips is accompanied by significant fattening in the 
Z-direction. 




Fig. 14. p = .1 : a classic dendrite. 




Fig. 15. p = .095 : fewer side branches. 
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Fig. 16. p = .09 : no significant side branches on this scale. 




Fig. 17. p = .082 : the tip undergoes a sandwich instabihty. 

The crystal in Fig. 17 is captured at about time 60000. The series of close-ups in Fig. 18 
provides another illustration of the sandwich instability — snapshots of the same snowfake are 
shown at time intervals of 1000, starting from time 37000. 
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Our final example, with p = .081, demonstrates that a further decrease in density makes the 
crystal increasingly three-dimensional. 



Fig. 19. Fattening from the tip inward at p = .081. 



9 Case study Hi : sandwich plates 

When growth in the Z-direction is much slower than in the xy-plane, outer ridges never develop. 
Instead, the dynamics grows a featureless prism, which, when sufficiently thick, undergoes a 
sandwich instability producing inner ridges. Much later the crystal experiences the branching 
instability, with plate-like branches that bear a superficial resemblance to Packard snowflakes 
|Pac[ IGG2j during early stages. 

Throughout the evolution the external surface of the crystal has few markings, whereas inside 
features include ridges and rifts, which signify gradual thinning of the plates from the center 
outward before the branching instability. 

The sole surface designs are reverse shapes, which occur when the crystal grows in the Z- 
direction from buds that arise close to the tips. These macrostep nuclei result in rapid growth of 
a single layer in the T-direction until this layer outlines a nearly circular hole near the crystal's 
center; the hole then proceeds to shrink much more slowly. 

We note that this observation provides a convincing explanation for the circular markings 
seen on many snow crystal photographs [Lib6l ILRj . It also suggests that ribs are predomi- 
nantly inner structures. While outer ribs could occur due to instabilities or changing conditions 
(cf. Fig. 11), there is scant evidence of them in electron microscope photographs [EMP], which 
completely obscure inner structure. On the other hand, those photos reveal an abundance of 
sandwich plates, which appear as the crystal centers, at the tips of the six main arms, and as 
side branches. 

We now present two examples. Both start from the canonical seed. In the first, depicted in 
Fig. 20, f3oi = 6, (3io = /32o = 2.5, = 2, and the remaining /3's are 1. All /^'s are .1, except 
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that Koi = .5, jii = .0001, and p = .08. The final radius of the crystal at the capture time 100000 
is about 150. Note that the main ridge is interrupted: while initially it connects the two plates 
(and it has darker color in the ray-traced image as the background can be seen through it), it 
later splits and each plate has its own ridge. There is a suggestion of this phenomenon in real 
crystals (e.g., on p. 26 of |Lib6| ). 




Fig. 20. A sandwich plate. 



Our second example (Figs. 21 and 22) has interrupted main ridges and a few ribs. The 
parameter set now has /3oi = 6.5, /3io = P20 = 2.7, and p = .15. The remaining values are as 
before, and the final sizes (this one at t = 36100) are comparable. We provide a few intermediate 
stages and a detail of the inner structure. Observe the buds at times 25883 and 31671; also note 
that the outermost rib at time 19000 later disappears. 




Fig. 21. Another sandwich plate. 
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Fig. 22. The plate of Fig. 21 at t = 19000,25883,25900,25950,26000,31671. The detail is 
from the first time, obtained by cutting the crystal along the plane z = and zooming in on 
the bottom half of the upper portion. 



10 Case study iv: the roles of drift and melting 

From some of the electron micros; raphs at |EMPj . it appears possible that the basal facets may 
have ridges and other markings on one side only, while the other side is nearly featureless. As 
far as we are aware, no attempt has been made to "turn over" these specimens and confirm 
the asymmetry, but |NK^ |Nel] offer a theoretical explanation. They suggest that the one-sided 
structure is a consequence of early growth and that ridges are actually vestiges of the skeleton 
of hollow prisms such as Fig. 31 in Section 11 (see Fig. 3 of [Nelj). In fact, it is widely held that 
the micron-scale prism from which a prototypical snowflake evolves develops slight asymmetries 
in the radii of its two basal facets, and that the larger facet acquires an increasing advantage 
from the feedback effect of diffusion-limited growth. As a result many crystals have a stunted 
hexagonal plate at their center. In |Nak| this effect is described on p. 206 and in sketch 15 of 
Fig. 369. 

Another potential source of asymmetry in the Z-direction is identified in Section 3.5 of |Iwaj 
and on p. 18 of jTEWFj, based on cloud tunnel experiments in the laboratory. Planar snowflakes 
evidently assume a preferred orientation parallel to the ground as they slowly fall, resulting in 
a small upward drift of the diffusion field relative to the crystal. 

We emulate these aspects of asymmetric growth by means of the drift (j) in step (Ic) of 
our algorithm and asymmetry of the initial seed as mentioned in Section 3. Consider first the 
snowfake of Fig. 1 and the closely related sectored plate in Fig. 23. The former starts from 
our fundamental prism and never undergoes the sandwich instability, but develops ridges on 
the bottom side and an almost featureless top due to the presence of = .01. The dynamic 
parameters of the sectored plate below are identical, but growth starts from a mesoscopic prism 
that is 5 cells high, with radius 7 at the top and 3 at the bottom. The idea here is to mimic 
the situation where the upper basal plate has established an advantage over the lower basal 
plate early in the evolution. As is clear from the side view, in contrast to Fig. 6, growth of 
the lower facet stops completely due to diffusion limitation even though the small drift offers 
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a slight advantage in the early stages. (According to |Iwaj , falling snowflakes prefer the more 
aerodynamically stable orientation of Fig. 23.) Very many photos of physical snow crystals 
show evidence of such a stunted simple plate at the center; see ^Lib6) . pp. 75-76, for further 
discussion. 




Fig. 23. A sectored plate with a stunted double, from the top (left) and side (right). 




Fig. 24. A fern dendrite for /iio = /i20 — -005. 

The remaining examples of this section also start from slightly asymmetric seeds, experience 
a small drift, and have almost all their external markings on one side. Our goal is to explore 
the role of the melting rate, in much the same way we studied density dependence in Section 
7, by varying /i in a series of snowfakes with all other parameters held fixed. In each instance, 
the seed has height 3, lower radius 2, and upper radius 1. For the next four crystals, /3oi = 3, 
/?io = P20 = Pu = 1.4, Pso = P21 = /331 = 1, = .1, = .01, and p = .14. Moreover /xqi = .002, 
M30 = Mil = M21 = M31 = .001 and we vary only the common value of pio = ^20. This value 
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governs the speed of tips and — as explained in Section 5 — has more effect in regions of low 
density, so an increase inhibits side branching. 

Like the sectored plates just discussed, these are relatively rare snowfakes with outside ridges 
on the main arms and most side branches. All our modeling experience suggests that crystal 
tips tend to symmetrize with respect to the T-direction, managing to avoid the sandwich in- 
stability only under quite special environmental conditions. We have seen little evidence in our 
simulations for the mechanism of ridge formation proposed in |NK[ INel ] , so we feel that drift 
is a more likely explanation of one-sided structures in snowflakes. 




Fig. 25. Reduced side branching for /j^iq — /X20 = .008. 




Fig. 26. Further reduction in the number of side branches for /xio = /i20 — -009. 
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Starting with the classic fern of Fig. 24, the common prism facet melting threshold /xio = 
is gradually increased to twice the original value in Figs. 25-7. Stellar dendrites with fewer and 
fewer side branches result, until the final snowfake has only a few short sandwich plates on the 
sides of each arm. 




Fig. 27. When /xio = /i20 — -01, very few side branches remain. 



The final example of this section is a classic simple star, a crystal with no side branches at all 
and a characteristic parabolic shape to its tips (cf. jLib6| . p. 57 bottom). This elegant snowfake 
required considerable tweaking of parameters; they are: /3oi = 3.1, /3io = 1.05, /320 — 1-03, 
Pii 1.04, Pso 1.02, (321 1.01, Psi = 1, = .01, /xoi M30 Mil M21 M3i .01 , 
Mio — M20 — .03, (j) — .005, and p — .16. 




Fig. 28. A simple star. 
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11 Case study v : needles and columns 

Let us now turn to the common but less familiar snow crystals that expand primarily in the Z- 
direction. As one would expect, these have /3oi small compared to f3io and f320, but surprisingly 
small advantage often suffices. We offer three snowfakes that emulate their physical counterparts 
quite well. All start from the canonical seed. Our first example, with a substantial bias toward 
attachment on the basal facets, is a (simple) needle. In Fig. 29, /3oi = 2, (3io = /320 = Pii = 4, 
/^ao — /?2i — /?3i = 1, = .1, = .001, = 0, and p = .1. This snowfake reproduces structure 
observed in nature and the laboratory: slender hollow tubes, often with cable-like protuberances 
at the ends (cf. Fig. 135 of [Nakj . pp. 67-68 of pb6] ). 




Fig. 29. A needle. 




Fig. 30. A hollow column. 
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Next, Fig. 30 simulates the common type of snow crystal known as a hollow column. Here 
the bias toward attachment on the basal facets is not as pronounced. The parameter set is: 
/?oi = 1, Pio = /320 = 2 Pso = Pii = (321 = .5 Psi = 1, = .1, /X = .01, = 0, and p = .1. 
Evidently, the hole starts developing early on. See pp. 64-66 of [Lib6j) for photos of actual 
hollow columns and a qualitative description of their growth. 

The final example of this section is a column whose facets are hollow as well. The morphology 
of Fig. 31 occurs when the rates of expansion in the two directions are not very different. Photos 
and a description of this sort of snowflake appear on pp. 35-37 of [Lib6j. Here /3oi = 1.5, 
Ao = /320 = 1-6 Pu = Pso = P21 = /331 = 1, ^ = -1, M = -015, = 0, and p = .1. 




12 Case study vi : change of environment 

In his pioneering research, Nakaya [NakJ reproduced several of the most striking types found in 
nature by subjecting the cold chamber in his lab to a precisely controlled schedule of temperature 
and humidity changes, either sudden or gradual. Based on such experiments, he argued that 
plates with dendritic extensions, for example, are formed when a snowflake's early growth occurs 
in the upper atmosphere and then it drops to another layer more conducive to branching ([NakJ, 
p. 16). 

In this section we mimic such varying environments by consider the effect of an abrupt change 
of parameters on some of our previous snowfakes. Let us begin with two examples of the type 
cited in the last paragraph: plates with dendritic extensions. Both start from a prism that is 
3 cells high with radius 2 at the top and 1 at the bottom. The first stage for both is a simple 
plate similar to the snowfake of Fig. 1, but with a delayed branching instability. The initial 
parameters are: /3oi 3.5, /3io /32o Ai 2.25, f^so (32i (^31 = 1, = .005, p = .001, 
= .01, and p = .12. The first stage runs until time 8000 in the first example, and until time 
12500 in the second. At that time most parameters remain the same, but in order to promote 
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branching we change Piq = /32o = Pii to 1.15 (resp. 1.4) and jj^iq = to .006 (resp. 004). 
The results, once the two snowfakes have reached a radius of 200 cells, are shown in Figs. 32-3. 
Predictably, the first example has more branching in its dendritic phase since the prism facet 
attachment threshold is lower. The large image on the cover of [LibGj shows a beautiful natural 
example of this type. 




Fig. 32. A plate with fern extensions. 




Fig. 33. A plate with dendrite extensions. 



A hybrid evolution at the opposite end of the spectrum is described in |Lib6| . pp. 51- 
53, and many of the most striking snowflakes in |LRj are of this type. As presumably in 
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nature, conditions need to be just right for the corresponding snowfake to evolve. In this 
vein, we present three snowfakes that begin as stellar dendrites with minimal branching and 
later encounter an environment promoting plates. All start from a prism of height 5 with top 
radius 6 and bottom radius 2. The first stage runs the simple star dynamics of Fig. 28 until 
time 4000, 3000, or 2000, respectively. Then new parameters for the three experiments with 
higher attachment thresholds are run until time, respectively, 24000, 20000, and close to 20000. 
Common parameters are: f3so — f^si = 1, = .1, p = .16. In Fig. 34, the remaining parameters 
are /3oi = 3.0, /?io = /320 = 2.2, (3u = 2.0, /32i = 1.1, fi = .01, (f) = .005. Note that in this 
instance the branches of the star broaden considerably after the change of environment, and the 
tips form sandwich plates. 




Fig. 34. A broad-branched stellar crystal with sandwich-plate extensions. 




Fig. 35. A broad-branched stellar crystal with sectored-plate extensions. 
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By raising the attachment thresholds somewhat we avoid the sandwich instabihty and obtain 
instead the sectored-plate extensions with outside ridges seen in Fig. 35. Here /3oi = 3.5, 
/3io = /320 = 2.45, (3ii = 2.25, /?2i = 1.1, m = m = -002, /x = .001 otherwise, = .015. 

Our final broad-branched example interpolates between the previous two. The values of (3 
are large enough to avoid the sandwich instability, but small enough that side branching leads 
to sectored plate structure of the extensions. Here /3oi = 3.0, Pio = P20 = 2.25, Pu = 2.05, 

= 1.05, /x = .001, = .015. 



Fig. 36. Another broad-branched stellar crystal. 

We conclude this case study with two crystals that combine a three-dimensional column 
and two-dimensional plates. These are the tsuzumi, or capped columns, described on pp. 69- 
74 of |Lib6 ] . They are thought to arise when crystals are transported to higher and colder 
regions of the atmosphere by a passing storm. Without a preferred orientation, it is most 
reasonable to model these as driftless. Both our snowfakes use the canonical seed and evolve 
with the parameters for the hollow column of Fig. 30 until time 20000. Then they run with 
new parameters that promote planar growth, until time 80000 for the first example, 60000 for 
the second. Common values for the two examples are: /3oi — 5, f3so — [^21 — f^si = 1, = .1, 
/i = .001, = 0, and p = .1. The difference is the common value /3io = P20 = Pu is 2.4 in 
Fig. 37, and 2.1 in Fig. 38. Higher attachment thresholds delay the branching instability in the 
first capped column so the caps are simple plates, as opposed to sectored plates in the second. 

The transition period from column to cap in lab tsuzumi is described in some detail by 
Nakaya ([Nak], p. 221; see also the sketch on p. 222). We remark that our snowfake versions 
evolve in the same way. Namely, for a considerable time after the change of environment, 
outward growth occurs almost exclusively along the 18 edges of the hexagonal column. This is 
a diffusion-limited effect similar to the hollowing in Fig. 31. Then, rather suddenly, growth in 
the T-direction takes over. 
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Fig. 37. A column capped with hexagonal plates. 




Fig. 38. A column capped with sectored plates. 



13 Case study vii : eccentric crystals 

This last section features snowfakes that result from a careful search through parameter space 
and are quite sensitive to any change. They are close-to-critical, near the phase boundary 
between dominant growth in the Z-direction and the T-direction. Consequently they may be 
rare in nature, though variants of some of the forms have been observed, and even represent 
morphological types in the Magono-Lee classification [MLJ . All our final examples start from 
the canonical seed. 
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As mentioned in Section 2, starting from a single cell our algorithm has a strong tendency 
to grow rapidly in the Z-direction due to the immediate onset of a needle instability. Even 
if the initial mesoscopic prism is wider in the T-direction, it is still quite common for this 
instability to arise later on if the dynamics are close to critical. After an initial phase of typical 
planar growth, needles suddenly nucleate at concentric locations scattered over the central plate 
or arms. Fig. 137 of |Nak| shows an excellent example of this type in nature, and our first 
two examples illustrate a similar phenomenon in our model. The conventional explanation for 
such hybrid types, called stellar crystals with needles in jMLj , involves a sudden change in the 
environment, but this is one of several cases where our algorithm suggests that homogeneous 
conditions can sometimes produce the same effect. 

Fig. 39 has features like a classic planar snowflake that has developed rime from attachment 
of surrounding water droplets. In fact these protrusions are potential needle instabilities — the 
two symmetric rings close to the center and the tips are stunted needles, whereas the intermediate 
needles have successfully nucleated. The parameters of this snowfake are: /3oi = 1.58, /3io = 
/?20 — Pii — 1-5, (3so — (321 — (^31 = 1, = .1, /X = .006, = and p — .1. Partial symmetry 
of bumps in many natural crystals, statistically unlikely to be the result of rime, often indicates 
vestiges of rims and ribs after sublimation, but can also be due to nascent needles, as in the 
middle specimen of Plate 116 in |Nak ] . Since the locations where needles nucleate are quite 
sensitive to changes in parameters, residual randomness in the mesoscopic dynamics is apt to 
degrade the symmetry. 





The next three examples have = 1, ^ = .03, a^io — 1^20 = -l? ^30 — -05, and /^n — hi2i — 
K31 = .01. The remaining parameters for Fig. 40 are kqi = .11 and p = .06. This snowfake is a 
rather extreme instance of a stellar crystal with needles in which the planar portion is a thick 
but very narrow simple star. 
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Our next two examples seem never to have been seen at all, and it is clear why: even if 
they managed to grow, their thin plates would be extremely brittle and susceptible to random 
fluctuations. They are characterized by very small differences in the growth rates. After starting 
as planar crystals, they suddenly nucleate thin structures extending into the third dimension. 
In Fig. 41 tvQi = .12 and p = .057; in Fig. 42 /^oi — -116 and p = .06. For obvious reasons, we 
call these hutter flakes. They are idealizations of the stellar crystals with spatial plates in |ML| : 
chaotic snow crystals with thin plates growing every which way are relatively common. 




Fig. 41. A butterflake with wings in the directions of the main arms. 
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Fig. 42. A butterflake with side wings. 



We conclude the paper with a family of five related examples. The first is a common sandwich 
plate (cf. p. 39, lower right, in [L ibGj ) with parameter values /3oi = 1.41, /?io = /?2o = 1-2 
Ai = /?30 = P2i = /331 = 1, = .1, /i = .025, = 0, and p = .09. 




Fig. 43. A sandwich plate with broad branches. 



The remaining four are minor perturbations, which nevertheless look quite different. Namely, 
even though their model parameters are constant over time, they undergo "exploding tips" quite 
similar to crystals such as the one in Fig. 35 that results from inhomogeneous environmental 
conditions. The principle behind all four variants is the same: eventually, the growing tip thick- 
ens and slows down considerably. Usually this happens close to the beginning of the evolution 
(as, in fact, occurred in the dynamics leading to Fig. 43), so the snowfake is unremarkable. But 
with some experimentation we find cases when the onset of the sandwich instability is delayed 
and the final picture can be quite dramatic. The complex inner patterns are the result of ex- 
traordinarily intricate dynamics. Parameter values that differ from those of Fig. 43 are given in 
the captions. 
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Fig. 46. Perturbed parameter: /3oi = 1.19. 




Fig. 47. Perturbed parameter: /3oi = 1.25. 
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